#clean up
rm(list=ls())
options(digits=8,scipen=8)

#Set Working Directory: REVISE THIS TO WHERE YOU SAVE THE DATA.
setwd('/Volumes/MONOGAN/scVice/data/monthly')

#Load Libraries
library(MASS)
library(TSA)
library(foreign)
library(xtable)
library(tseries)
library(vars)
library(car)
library(XLConnect)

#load files
#comprehensive<-read.csv('comprehensiveData8109.csv')
comprehensive<-readWorksheetFromFile('comprehensive_data_8.10.16.xls',sheet=1,header=TRUE)

##CREATE PRESIDENTIAL AND CONGRESSIONAL VARIABLES##
comprehensive$reagan<-as.numeric(comprehensive$year>=1981 & comprehensive$year<=1988)
comprehensive$bush41<-as.numeric(comprehensive$year>=1989 & comprehensive$year<=1992)
comprehensive$clinton<-as.numeric(comprehensive$year>=1993 & comprehensive$year<=2000)
comprehensive$bush43<-as.numeric(comprehensive$year>=2001 & comprehensive$year<=2008)
comprehensive$obama<-as.numeric(comprehensive$year>=2009 & comprehensive$year<=2010)

comprehensive$tobaccoCONG<-comprehensive$c_bills_tobacco+comprehensive$c_hearings_tobacco+comprehensive$c_laws_tobacco
comprehensive$tobaccoPRES<-comprehensive$p_pressconf_tobacco+comprehensive$p_execord_tobacco+comprehensive$p_radio_tobacco+comprehensive$tobacco_presTV

comprehensive$alcoholCONG<-comprehensive$c_bills_alcohol+comprehensive$c_hearings_alcohol+comprehensive$c_laws_alcohol
comprehensive$alcoholPRES<-comprehensive$p_pressconf_alcohol+comprehensive$p_execord_alcohol+comprehensive$p_radio_alcohol+comprehensive$alcohol_presTV

comprehensive$p_pressconf_porn<-0
comprehensive$pornCONG<-comprehensive$c_bills_porn+comprehensive$c_hearings_porn+comprehensive$c_laws_porn
comprehensive$pornPRES<-comprehensive$p_pressconf_porn+comprehensive$p_execord_porn+comprehensive$p_radio_porn+comprehensive$porn_presTV

comprehensive$gamblingCONG<-comprehensive$c_bills_gambling+comprehensive$c_hearings_gambling+comprehensive$c_laws_gambling
comprehensive$gamblingPRES<-comprehensive$p_pressconf_gambling+comprehensive$p_execord_gambling+comprehensive$p_radio_gambling+comprehensive$gambling_presTV

comprehensive$drugsCONG<-comprehensive$c_bills_drugs+comprehensive$c_hearings_drugs+comprehensive$c_laws_drugs
comprehensive$drugsPRES<-comprehensive$p_pressconf_drugs+comprehensive$p_execord_drugs+comprehensive$p_radio_drugs+comprehensive$illegal_drugs_presTV

comprehensive$cocaineCONG<-comprehensive$c_bills_cocaine+comprehensive$c_hearings_cocaine+comprehensive$c_laws_cocaine
comprehensive$cocainePRES<-comprehensive$p_pressconf_cocaine+comprehensive$p_execord_cocaine+comprehensive$p_radio_cocaine+comprehensive$cocaine_presTV

comprehensive$c_laws_marj<-0
comprehensive$marjCONG<-comprehensive$c_bills_marj+comprehensive$c_hearings_marj+comprehensive$c_laws_marj
comprehensive$marjPRES<-comprehensive$p_pressconf_marj+comprehensive$p_execord_marj+comprehensive$p_radio_marj+comprehensive$marijuana_presTV




####DESCRIPTIVE STATISTICS#####
names(comprehensive)
desc<-subset(comprehensive,select=c(mood,m_tobacco,tobaccoCONG,tobaccoPRES,m_drugs,drugsCONG,drugsPRES))
xtable(cbind(apply(desc,2,mean),apply(desc,2,sd),apply(desc,2,min),apply(desc,2,max)),digits=4)

##GRAPHS##
#pdf('../instAttention.pdf')
attention<-matrix(NA,nrow=7,ncol=2)
attention[1,1]<-sum(comprehensive$tobaccoCONG)
attention[1,2]<-sum(comprehensive$tobaccoPRES)
attention[2,1]<-sum(comprehensive$alcoholCONG)
attention[2,2]<-sum(comprehensive$alcoholPRES)
attention[3,1]<-sum(comprehensive$pornCONG)
attention[3,2]<-sum(comprehensive$pornPRES)
attention[4,1]<-sum(comprehensive$gamblingCONG)
attention[4,2]<-sum(comprehensive$gamblingPRES)
attention[5,1]<-sum(comprehensive$drugsCONG)
attention[5,2]<-sum(comprehensive$drugsPRES)
attention[6,1]<-sum(comprehensive$cocaineCONG)
attention[6,2]<-sum(comprehensive$cocainePRES)
attention[7,1]<-sum(comprehensive$marjCONG)
attention[7,2]<-sum(comprehensive$marjPRES)
rownames(attention)<-c('Tobacco','Alcohol','Pornography','Gambling','Drugs','Cocaine','Marijuana')
colnames(attention)<-c('Congress','President')
attention<-attention[order(attention[,1]),]
par(mar=c(5.1,10,4.1,2.1))
barplot(t(attention),beside=T,horiz=T,las=1,xlab='Number of Actions',xlim=c(0,2150),legend=c('Congress','President'),args.legend=list(x='bottomright'))
box()
#dev.off()

#pdf('mediaCoverage.pdf')
plot(comprehensive$m_tobacco,type='l',axes=F,xlab="Month",ylab="Number of Stories",ylim=c(0,300))
lines(x=1:360,y=comprehensive$m_drugs,lty=2,col='blue')
axis(2)
axis(1,at=c(1,61,121,181,241,301),labels=c(1981,1986,1991,1996,2001,2006))
box()
text(x=150,y=250,"Drugs")
text(x=150,y=10,"Tobacco")
arrows(x0=145,x1=115,y0=245,y1=225,length=.1)
arrows(x0=145,x1=122,y0=15,y1=28,length=.1)
#dev.off()




####TOBACCO MODEL#####
#create subsets
tobacco<-subset(comprehensive, select=c(m_tobacco,tobaccoCONG,tobaccoPRES))

#estimate VAR 
tobacco.VAR<-VAR(tobacco, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$sc_fda_tobacco,comprehensive$sc_lorillard_tobacco)); tobacco.VAR; AIC(tobacco.VAR)

#assess the model
#plot(tobacco.VAR)
plot(irf(tobacco.VAR,boot=FALSE,impulse="m_tobacco",response="m_tobacco"))
plot(irf(tobacco.VAR,boot=FALSE,impulse="m_tobacco",response=c("tobaccoCONG","tobaccoPRES")))
plot(irf(tobacco.VAR,boot=FALSE,impulse="tobaccoCONG",response="tobaccoCONG"))
plot(irf(tobacco.VAR,boot=FALSE,impulse="tobaccoCONG",response=c("m_tobacco","tobaccoPRES")))
plot(irf(tobacco.VAR,boot=FALSE,impulse="tobaccoPRES",response="tobaccoPRES"))
plot(irf(tobacco.VAR,boot=FALSE,impulse="tobaccoPRES",response=c("tobaccoCONG","m_tobacco")))

###specific tests###
#cases on outcomes
summary(tobacco.VAR$varresult$m_tobacco) 
summary(tobacco.VAR$varresult$tobaccoCONG) 
summary(tobacco.VAR$varresult$tobaccoPRES) 

#Items on Press Coverage
a<-linearHypothesis(tobacco.VAR$varresult$m_tobacco,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
b<-linearHypothesis(tobacco.VAR$varresult$m_tobacco,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
c<-linearHypothesis(tobacco.VAR$varresult$m_tobacco,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(tobacco.VAR$varresult$tobaccoCONG,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
e<-linearHypothesis(tobacco.VAR$varresult$tobaccoCONG,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
f<-linearHypothesis(tobacco.VAR$varresult$tobaccoCONG,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Items on President
g<-linearHypothesis(tobacco.VAR$varresult$tobaccoPRES,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
h<-linearHypothesis(tobacco.VAR$varresult$tobaccoPRES,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
i<-linearHypothesis(tobacco.VAR$varresult$tobaccoPRES,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Table of Granger tests
tobacco.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(tobacco.granger)<-colnames(tobacco.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(tobacco.granger,caption="\\textbf{Granger Causality Tests for Tobacco VAR Model}",label="tobacco.granger",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i

###Draw the Effect of FDA v. Brown & Williamson Tobacco Corp.###
##Mean Effect Only##
#Intervention: March 2000, month #231
#Values of endogenous variables the month before the intervention, February 2000
feb.00<-c(95,2,0)

#Define inputs and make predictions for two prior months and intervention month.
x.00<-c(rep(feb.00,12),1,229,62.312,0,0)
m.00<-x.00%*%tobacco.VAR$varresult$m_tobacco$coefficients
c.00<-x.00%*%tobacco.VAR$varresult$tobaccoCONG$coefficients
p.00<-x.00%*%tobacco.VAR$varresult$tobaccoPRES$coefficients

x.0<-c(rep(feb.00,12),1,230,62.312,0,0)
m.0<-x.0%*%tobacco.VAR$varresult$m_tobacco$coefficients
c.0<-x.0%*%tobacco.VAR$varresult$tobaccoCONG$coefficients
p.0<-x.0%*%tobacco.VAR$varresult$tobaccoPRES$coefficients

all.inp<-x.1<-c(rep(feb.00,12),1,231,62.312,1,0)
m.1<-x.1%*%tobacco.VAR$varresult$m_tobacco$coefficients
c.1<-x.1%*%tobacco.VAR$varresult$tobaccoCONG$coefficients
p.1<-x.1%*%tobacco.VAR$varresult$tobaccoPRES$coefficients
new.pred<-c(m.1,c.1,p.1)

#Use a loop to make predictions in subsequent months.
for(i in 2:10){
	all.inp<-assign(paste("x",i,sep="."),c(new.pred,all.inp[-c(34:41)],1,230+i,62.312,0,0))
	assign(paste("m",i,sep="."),get(paste("x",i,sep="."))%*%tobacco.VAR$varresult$m_tobacco$coefficients)
	assign(paste("c",i,sep="."),get(paste("x",i,sep="."))%*%tobacco.VAR$varresult$tobaccoCONG$coefficients)
	assign(paste("p",i,sep="."),get(paste("x",i,sep="."))%*%tobacco.VAR$varresult$tobaccoPRES$coefficients)
	new.pred<-c(get(paste("m",i,sep=".")),get(paste("c",i,sep=".")),get(paste("p",i,sep=".")))
}

#Save Mean Figures and get preliminary plots
Time<-c(-1:10)
t.m.forecast<-c(m.00,m.0,m.1,m.2,m.3,m.4,m.5,m.6,m.7,m.8,m.9,m.10)
plot(y=t.m.forecast,x=Time,ylab="Media Forecast",type='o')
t.c.forecast<-c(c.00,c.0,c.1,c.2,c.3,c.4,c.5,c.6,c.7,c.8,c.9,c.10)
plot(y=t.c.forecast,x=Time,ylab="Congressional Forecast",type='o')
t.p.forecast<-c(p.00,p.0,p.1,p.2,p.3,p.4,p.5,p.6,p.7,p.8,p.9,p.10)
plot(y=t.p.forecast,x=Time,ylab="Presidential Forecast",type='o')

##Create Clarify-Like Confidence Intervals##
t.m.clarify<-matrix(NA,nrow=1000,ncol=12)
t.c.clarify<-matrix(NA,nrow=1000,ncol=12)
t.p.clarify<-matrix(NA,nrow=1000,ncol=12)

#Clarify simulations
tobacco.media.sim<-mvrnorm(n=1000,mu=tobacco.VAR$varresult$m_tobacco$coefficients,Sigma=vcov(tobacco.VAR$varresult$m_tobacco),empirical=TRUE)
tobacco.cong.sim<-mvrnorm(n=1000,mu=tobacco.VAR$varresult$tobaccoCONG$coefficients,Sigma=vcov(tobacco.VAR$varresult$tobaccoCONG),empirical=TRUE)
tobacco.pres.sim<-mvrnorm(n=1000,mu=tobacco.VAR$varresult$tobaccoPRES$coefficients,Sigma=vcov(tobacco.VAR$varresult$tobaccoPRES),empirical=TRUE)

#Loop for simulated forecasts. Relies on standing terms for pre-intervention endogenous variables, x.00, x.0, and x.1.
for(j in 1:1000){
	tobacco.media.coef<-tobacco.media.sim[j,]
	tobacco.cong.coef<-tobacco.cong.sim[j,]
	tobacco.pres.coef<-tobacco.pres.sim[j,]

	t.m.clarify[j,1]<-m.00<-x.00%*%tobacco.media.coef
	t.c.clarify[j,1]<-c.00<-x.00%*%tobacco.cong.coef
	t.p.clarify[j,1]<-p.00<-x.00%*%tobacco.pres.coef

	t.m.clarify[j,2]<-m.0<-x.0%*%tobacco.media.coef
	t.c.clarify[j,2]<-c.0<-x.0%*%tobacco.cong.coef
	t.p.clarify[j,2]<-p.0<-x.0%*%tobacco.pres.coef

	all.inp<-x.1
	t.m.clarify[j,3]<-m.1<-x.1%*%tobacco.media.coef
	t.c.clarify[j,3]<-c.1<-x.1%*%tobacco.cong.coef
	t.p.clarify[j,3]<-p.1<-x.1%*%tobacco.pres.coef
	new.pred<-c(m.1,c.1,p.1)

	for(i in 2:10){
		all.inp<-assign(paste("x",i,sep="."),c(new.pred,all.inp[-c(34:41)],1,230+i,62.312,0,0))
		t.m.clarify[j,i+2]<-assign(paste("m",i,sep="."),get(paste("x",i,sep="."))%*%tobacco.media.coef)
		t.c.clarify[j,i+2]<-assign(paste("c",i,sep="."),get(paste("x",i,sep="."))%*%tobacco.cong.coef)
		t.p.clarify[j,i+2]<-assign(paste("p",i,sep="."),get(paste("x",i,sep="."))%*%tobacco.pres.coef)
		new.pred<-c(get(paste("m",i,sep=".")),get(paste("c",i,sep=".")),get(paste("p",i,sep=".")))
		}
}

#Lower and Upper Bounds
t.m.5<-apply(t.m.clarify,2, quantile, .05)
t.c.5<-apply(t.c.clarify,2, quantile, .05)
t.p.5<-apply(t.p.clarify,2, quantile, .05)
t.m.95<-apply(t.m.clarify,2, quantile, .95)
t.c.95<-apply(t.c.clarify,2, quantile, .95)
t.p.95<-apply(t.p.clarify,2, quantile, .95)

#Draw Figures with Confidence intervals
Time<-c(-1:10)

##pdf('../tobaccoMediaForecast.pdf')
plot(y=t.m.forecast,x=Time,ylab="Media Forecast",type='l',lwd=2,ylim=c(65,155),axes=F)
lines(y=t.m.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=t.m.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Feb 00','Apr 00','June 00','Aug 00','Sept 00','Dec 00'))
box()
##dev.off()

##pdf('../tobaccoCongForecast.pdf')
plot(y=t.c.forecast,x=Time,ylab="Congressional Forecast",type='l',lwd=2,ylim=c(-1,8),axes=F)
lines(y=t.c.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=t.c.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Feb 00','Apr 00','June 00','Aug 00','Sept 00','Dec 00'))
box()
##dev.off()

##pdf('../tobaccoPresForecast.pdf')
plot(y=t.p.forecast,x=Time,ylab="Presidential Forecast",type='l',lwd=2,ylim=c(-1,3),axes=F)
lines(y=t.p.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=t.p.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Feb 00','Apr 00','June 00','Aug 00','Sept 00','Dec 00'))
box()
##dev.off()




####DRUGS MODEL#####
#create subsets
drugs<-subset(comprehensive, select=c(m_drugs,drugsCONG,drugsPRES))

#estimate VAR 
drugs.VAR<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$sc_leon_drugs,comprehensive$sc_raab_drugs,comprehensive$sc_skinner_drugs,comprehensive$sc_smith_drugs,comprehensive$sc_austin_drugs,comprehensive$sc_vernonia_drugs,comprehensive$sc_ferguson_drugs,comprehensive$sc_kyllo_drugs,comprehensive$sc_earls_drugs,comprehensive$sc_benitez_drugs,comprehensive$sc_raich_drugs,comprehensive$sc_georgia_drugs,comprehensive$sc_hudson_drugs,comprehensive$sc_diaz_drugs)); drugs.VAR; AIC(drugs.VAR)

#assess the model
#plot(drugs.VAR)
plot(irf(drugs.VAR,boot=FALSE,impulse="m_drugs",response="m_drugs"))
plot(irf(drugs.VAR,boot=FALSE,impulse="m_drugs",response=c("drugsCONG","drugsPRES")))
plot(irf(drugs.VAR,boot=FALSE,impulse="drugsCONG",response="drugsCONG"))
plot(irf(drugs.VAR,boot=FALSE,impulse="drugsCONG",response=c("m_drugs","drugsPRES")))
plot(irf(drugs.VAR,boot=FALSE,impulse="drugsPRES",response="drugsPRES"))
plot(irf(drugs.VAR,boot=FALSE,impulse="drugsPRES",response=c("drugsCONG","m_drugs")))

###specific tests###
#cases on outcomes
summary(drugs.VAR$varresult$m_drugs) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
summary(drugs.VAR$varresult$drugsCONG) #Employment Division v. Smith (April 1990).
summary(drugs.VAR$varresult$drugsPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).
xtable(drugs.VAR$varresult$m_drugs) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
xtable(drugs.VAR$varresult$drugsCONG) #Employment Division v. Smith (April 1990).
xtable(drugs.VAR$varresult$drugsPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).

#Items on Press Coverage
a<-linearHypothesis(drugs.VAR$varresult$m_drugs,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
b<-linearHypothesis(drugs.VAR$varresult$m_drugs,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
c<-linearHypothesis(drugs.VAR$varresult$m_drugs,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(drugs.VAR$varresult$drugsCONG,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
e<-linearHypothesis(drugs.VAR$varresult$drugsCONG,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
f<-linearHypothesis(drugs.VAR$varresult$drugsCONG,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on President
g<-linearHypothesis(drugs.VAR$varresult$drugsPRES,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
h<-linearHypothesis(drugs.VAR$varresult$drugsPRES,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
i<-linearHypothesis(drugs.VAR$varresult$drugsPRES,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Table of Granger tests
drugs.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(drugs.granger)<-colnames(drugs.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(drugs.granger,caption="\\textbf{Granger Causality Tests for Drugs VAR Model}",label="drugs.granger",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i

###Draw the Effect of National Treasury Employees' Union v. Von Raab###
##Mean Effect Only##
#Intervention: March 1989, month #99
#Values of endogenous variables the month before the intervention, Feb 1989
feb.89<-c(132,6,1)

#Define inputs and make predictions for two prior months and intervention month.
x.00<-c(rep(feb.89,12),1,97,68.267,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
m.00<-x.00%*%drugs.VAR$varresult$m_drugs$coefficients
c.00<-x.00%*%drugs.VAR$varresult$drugsCONG$coefficients
p.00<-x.00%*%drugs.VAR$varresult$drugsPRES$coefficients

x.0<-c(rep(feb.89,12),1,98,68.267,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
m.0<-x.0%*%drugs.VAR$varresult$m_drugs$coefficients
c.0<-x.0%*%drugs.VAR$varresult$drugsCONG$coefficients
p.0<-x.0%*%drugs.VAR$varresult$drugsPRES$coefficients

all.inp<-x.1<-c(rep(feb.89,12),1,99,68.267,0,1,0,0,0,0,0,0,0,0,0,0,0,0)
m.1<-x.1%*%drugs.VAR$varresult$m_drugs$coefficients
c.1<-x.1%*%drugs.VAR$varresult$drugsCONG$coefficients
p.1<-x.1%*%drugs.VAR$varresult$drugsPRES$coefficients
new.pred<-c(m.1,c.1,p.1)

#Use a loop to make predictions in subsequent months.
for(i in 2:10){
	all.inp<-assign(paste("x",i,sep="."),c(new.pred,all.inp[-c(34:53)],1,98+i,68.267,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
	assign(paste("m",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$m_drugs$coefficients)
	assign(paste("c",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$drugsCONG$coefficients)
	assign(paste("p",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$drugsPRES$coefficients)
	new.pred<-c(get(paste("m",i,sep=".")),get(paste("c",i,sep=".")),get(paste("p",i,sep=".")))
}

#Save Mean Figures and get preliminary plots
Time<-c(-1:10)
d.m.forecast<-c(m.00,m.0,m.1,m.2,m.3,m.4,m.5,m.6,m.7,m.8,m.9,m.10)
plot(y=d.m.forecast,x=Time,ylab="Media Forecast",type='o')
d.c.forecast<-c(c.00,c.0,c.1,c.2,c.3,c.4,c.5,c.6,c.7,c.8,c.9,c.10)
plot(y=d.c.forecast,x=Time,ylab="Congressional Forecast",type='o')
d.p.forecast<-c(p.00,p.0,p.1,p.2,p.3,p.4,p.5,p.6,p.7,p.8,p.9,p.10)
plot(y=d.p.forecast,x=Time,ylab="Presidential Forecast",type='o')

##Create Clarify-Like Confidence Intervals##
d.m.clarify<-matrix(NA,nrow=1000,ncol=12)
d.c.clarify<-matrix(NA,nrow=1000,ncol=12)
d.p.clarify<-matrix(NA,nrow=1000,ncol=12)

#Clarify simulations
drugs.media.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$m_drugs$coefficients,Sigma=vcov(drugs.VAR$varresult$m_drugs),empirical=TRUE)
drugs.cong.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$drugsCONG$coefficients,Sigma=vcov(drugs.VAR$varresult$drugsCONG),empirical=TRUE)
drugs.pres.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$drugsPRES$coefficients,Sigma=vcov(drugs.VAR$varresult$drugsPRES),empirical=TRUE)

#Loop for simulated forecasts. This calls on standing terms defined in the mean forecast: Start endogenous variables, x00, x0, and x1.
for(j in 1:1000){
	drugs.media.coef<-drugs.media.sim[j,]
	drugs.cong.coef<-drugs.cong.sim[j,]
	drugs.pres.coef<-drugs.pres.sim[j,]

	d.m.clarify[j,1]<-m.00<-x.00%*%drugs.media.coef
	d.c.clarify[j,1]<-c.00<-x.00%*%drugs.cong.coef
	d.p.clarify[j,1]<-p.00<-x.00%*%drugs.pres.coef

	d.m.clarify[j,2]<-m.0<-x.0%*%drugs.media.coef
	d.c.clarify[j,2]<-c.0<-x.0%*%drugs.cong.coef
	d.p.clarify[j,2]<-p.0<-x.0%*%drugs.pres.coef

	all.inp<-x.1
	d.m.clarify[j,3]<-m.1<-x.1%*%drugs.media.coef
	d.c.clarify[j,3]<-c.1<-x.1%*%drugs.cong.coef
	d.p.clarify[j,3]<-p.1<-x.1%*%drugs.pres.coef
	new.pred<-c(m.1,c.1,p.1)

	for(i in 2:10){
		all.inp<-assign(paste("x",i,sep="."),c(new.pred,all.inp[-c(34:53)],1,98+i,68.267,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
		d.m.clarify[j,i+2]<-assign(paste("m",i,sep="."),get(paste("x",i,sep="."))%*%drugs.media.coef)
		d.c.clarify[j,i+2]<-assign(paste("c",i,sep="."),get(paste("x",i,sep="."))%*%drugs.cong.coef)
		d.p.clarify[j,i+2]<-assign(paste("p",i,sep="."),get(paste("x",i,sep="."))%*%drugs.pres.coef)
		new.pred<-c(get(paste("m",i,sep=".")),get(paste("c",i,sep=".")),get(paste("p",i,sep=".")))
		}
}

#Lower and Upper Bounds
d.m.5<-apply(d.m.clarify,2, quantile, .05)
d.c.5<-apply(d.c.clarify,2, quantile, .05)
d.p.5<-apply(d.p.clarify,2, quantile, .05)
d.m.95<-apply(d.m.clarify,2, quantile, .95)
d.c.95<-apply(d.c.clarify,2, quantile, .95)
d.p.95<-apply(d.p.clarify,2, quantile, .95)

#Draw Figures with Confidence intervals
Time<-c(-1:10)

##pdf('../drugsMediaForecast.pdf')
plot(y=d.m.forecast,x=Time,ylab="Media Forecast",type='l',lwd=2,ylim=c(122,262),axes=F)
lines(y=d.m.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.m.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Feb 89','Apr 89','June 89','Aug 89','Oct 89','Dec 89'))
box()
##dev.off()

##pdf('../drugsCongForecast.pdf')
plot(y=d.c.forecast,x=Time,ylab="Congressional Forecast",type='l',lwd=2,ylim=c(-8.5,20),axes=F)
lines(y=d.c.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.c.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Feb 89','Apr 89','June 89','Aug 89','Oct 89','Dec 89'))
box()
##dev.off()

##pdf('../drugsPresForecast.pdf')
plot(y=d.p.forecast,x=Time,ylab="Presidential Forecast",type='l',lwd=2,ylim=c(0,4.5),axes=F)
lines(y=d.p.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.p.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Feb 89','Apr 89','June 89','Aug 89','Oct 89','Dec 89'))
box()
##dev.off()


###Draw the Effect of Employment Division v. Smith (1990)###
##Mean Effect Only##
#Intervention: April 1990, month #112
#Values of endogenous variables the month before the intervention, March 1990
#comprehensive[111:113,1:10]
mar.90<-c(160,19,2)

#Define inputs and make predictions for two prior months and intervention month.
x.00<-c(rep(mar.90,12),1,110,64.982,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
m.00<-x.00%*%drugs.VAR$varresult$m_drugs$coefficients
c.00<-x.00%*%drugs.VAR$varresult$drugsCONG$coefficients
p.00<-x.00%*%drugs.VAR$varresult$drugsPRES$coefficients

x.0<-c(rep(mar.90,12),1,111,64.982,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
m.0<-x.0%*%drugs.VAR$varresult$m_drugs$coefficients
c.0<-x.0%*%drugs.VAR$varresult$drugsCONG$coefficients
p.0<-x.0%*%drugs.VAR$varresult$drugsPRES$coefficients

all.inp<-x.1<-c(rep(mar.90,12),1,112,64.982,0,0,0,1,0,0,0,0,0,0,0,0,0,0)
m.1<-x.1%*%drugs.VAR$varresult$m_drugs$coefficients
c.1<-x.1%*%drugs.VAR$varresult$drugsCONG$coefficients
p.1<-x.1%*%drugs.VAR$varresult$drugsPRES$coefficients
new.pred<-c(m.1,c.1,p.1)

#Use a loop to make predictions in subsequent months.
for(i in 2:10){
	all.inp<-assign(paste("x",i,sep="."),c(new.pred,all.inp[-c(34:53)],1,111+i,64.982,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
	assign(paste("m",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$m_drugs$coefficients)
	assign(paste("c",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$drugsCONG$coefficients)
	assign(paste("p",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$drugsPRES$coefficients)
	new.pred<-c(get(paste("m",i,sep=".")),get(paste("c",i,sep=".")),get(paste("p",i,sep=".")))
}

#Save Mean Figures and get preliminary plots
Time<-c(-1:10)
d.m.forecast<-c(m.00,m.0,m.1,m.2,m.3,m.4,m.5,m.6,m.7,m.8,m.9,m.10)
plot(y=d.m.forecast,x=Time,ylab="Media Forecast",type='o')
d.c.forecast<-c(c.00,c.0,c.1,c.2,c.3,c.4,c.5,c.6,c.7,c.8,c.9,c.10)
plot(y=d.c.forecast,x=Time,ylab="Congressional Forecast",type='o')
d.p.forecast<-c(p.00,p.0,p.1,p.2,p.3,p.4,p.5,p.6,p.7,p.8,p.9,p.10)
plot(y=d.p.forecast,x=Time,ylab="Presidential Forecast",type='o')

##Create Clarify-Like Confidence Intervals##
d.m.clarify<-matrix(NA,nrow=1000,ncol=12)
d.c.clarify<-matrix(NA,nrow=1000,ncol=12)
d.p.clarify<-matrix(NA,nrow=1000,ncol=12)

#Clarify simulations
drugs.media.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$m_drugs$coefficients,Sigma=vcov(drugs.VAR$varresult$m_drugs),empirical=TRUE)
drugs.cong.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$drugsCONG$coefficients,Sigma=vcov(drugs.VAR$varresult$drugsCONG),empirical=TRUE)
drugs.pres.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$drugsPRES$coefficients,Sigma=vcov(drugs.VAR$varresult$drugsPRES),empirical=TRUE)

#Loop for simulated forecasts. This calls on standing terms defined in the mean forecast: Start endogenous variables, x00, x0, and x1.
for(j in 1:1000){
	drugs.media.coef<-drugs.media.sim[j,]
	drugs.cong.coef<-drugs.cong.sim[j,]
	drugs.pres.coef<-drugs.pres.sim[j,]

	d.m.clarify[j,1]<-m.00<-x.00%*%drugs.media.coef
	d.c.clarify[j,1]<-c.00<-x.00%*%drugs.cong.coef
	d.p.clarify[j,1]<-p.00<-x.00%*%drugs.pres.coef

	d.m.clarify[j,2]<-m.0<-x.0%*%drugs.media.coef
	d.c.clarify[j,2]<-c.0<-x.0%*%drugs.cong.coef
	d.p.clarify[j,2]<-p.0<-x.0%*%drugs.pres.coef

	all.inp<-x.1
	d.m.clarify[j,3]<-m.1<-x.1%*%drugs.media.coef
	d.c.clarify[j,3]<-c.1<-x.1%*%drugs.cong.coef
	d.p.clarify[j,3]<-p.1<-x.1%*%drugs.pres.coef
	new.pred<-c(m.1,c.1,p.1)

	for(i in 2:10){
		all.inp<-assign(paste("x",i,sep="."),c(new.pred,all.inp[-c(34:53)],1,111+i,64.982,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
		d.m.clarify[j,i+2]<-assign(paste("m",i,sep="."),get(paste("x",i,sep="."))%*%drugs.media.coef)
		d.c.clarify[j,i+2]<-assign(paste("c",i,sep="."),get(paste("x",i,sep="."))%*%drugs.cong.coef)
		d.p.clarify[j,i+2]<-assign(paste("p",i,sep="."),get(paste("x",i,sep="."))%*%drugs.pres.coef)
		new.pred<-c(get(paste("m",i,sep=".")),get(paste("c",i,sep=".")),get(paste("p",i,sep=".")))
		}
}

#Lower and Upper Bounds
d.m.5<-apply(d.m.clarify,2, quantile, .05)
d.c.5<-apply(d.c.clarify,2, quantile, .05)
d.p.5<-apply(d.p.clarify,2, quantile, .05)
d.m.95<-apply(d.m.clarify,2, quantile, .95)
d.c.95<-apply(d.c.clarify,2, quantile, .95)
d.p.95<-apply(d.p.clarify,2, quantile, .95)

#Draw Figures with Confidence intervals
Time<-c(-1:10)

##pdf('drugsMediaForecastB.pdf')
plot(y=d.m.forecast,x=Time,ylab="Media Forecast",type='l',lwd=2,ylim=c(min(d.m.5),max(d.m.95)),axes=F)
lines(y=d.m.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.m.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Mar 90','May 90','July 90','Sept 90','Nov 90','Jan 91'))
box()
##dev.off()
min(d.m.5);max(d.m.95)

##pdf('drugsCongForecastB.pdf')
plot(y=d.c.forecast,x=Time,ylab="Congressional Forecast",type='l',lwd=2,ylim=c(min(d.c.5),max(d.c.95)),axes=F)
lines(y=d.c.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.c.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Mar 90','May 90','July 90','Sept 90','Nov 90','Jan 91'))
box()
##dev.off()

##pdf('drugsPresForecastB.pdf')
plot(y=d.p.forecast,x=Time,ylab="Presidential Forecast",type='l',lwd=2,ylim=c(min(d.p.5),max(d.p.95)),axes=F)
lines(y=d.p.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.p.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Mar 90','May 90','July 90','Sept 90','Nov 90','Jan 91'))
box()
##dev.off()


###Draw the Effect of Vernonia School District v. Acton (1995)###
##Mean Effect Only##
#Intervention: June 1995, month #174
#Values of endogenous variables the month before the intervention, May 1995
#comprehensive[173:174,1:10]
may.95<-c(99,2,1)

#Define inputs and make predictions for two prior months and intervention month.
x.00<-c(rep(may.95,12),1,172,55.398,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
m.00<-x.00%*%drugs.VAR$varresult$m_drugs$coefficients
c.00<-x.00%*%drugs.VAR$varresult$drugsCONG$coefficients
p.00<-x.00%*%drugs.VAR$varresult$drugsPRES$coefficients

x.0<-c(rep(may.95,12),1,173,55.398,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
m.0<-x.0%*%drugs.VAR$varresult$m_drugs$coefficients
c.0<-x.0%*%drugs.VAR$varresult$drugsCONG$coefficients
p.0<-x.0%*%drugs.VAR$varresult$drugsPRES$coefficients

all.inp<-x.1<-c(rep(may.95,12),1,174,55.398,0,0,0,0,0,1,0,0,0,0,0,0,0,0)
m.1<-x.1%*%drugs.VAR$varresult$m_drugs$coefficients
c.1<-x.1%*%drugs.VAR$varresult$drugsCONG$coefficients
p.1<-x.1%*%drugs.VAR$varresult$drugsPRES$coefficients
new.pred<-c(m.1,c.1,p.1)

#Use a loop to make predictions in subsequent months.
for(i in 2:10){
	all.inp<-assign(paste("x",i,sep="."),c(new.pred,all.inp[-c(34:53)],1,173+i,55.398,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
	assign(paste("m",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$m_drugs$coefficients)
	assign(paste("c",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$drugsCONG$coefficients)
	assign(paste("p",i,sep="."),get(paste("x",i,sep="."))%*%drugs.VAR$varresult$drugsPRES$coefficients)
	new.pred<-c(get(paste("m",i,sep=".")),get(paste("c",i,sep=".")),get(paste("p",i,sep=".")))
}

#Save Mean Figures and get preliminary plots
Time<-c(-1:10)
d.m.forecast<-c(m.00,m.0,m.1,m.2,m.3,m.4,m.5,m.6,m.7,m.8,m.9,m.10)
plot(y=d.m.forecast,x=Time,ylab="Media Forecast",type='o')
d.c.forecast<-c(c.00,c.0,c.1,c.2,c.3,c.4,c.5,c.6,c.7,c.8,c.9,c.10)
plot(y=d.c.forecast,x=Time,ylab="Congressional Forecast",type='o')
d.p.forecast<-c(p.00,p.0,p.1,p.2,p.3,p.4,p.5,p.6,p.7,p.8,p.9,p.10)
plot(y=d.p.forecast,x=Time,ylab="Presidential Forecast",type='o')

##Create Clarify-Like Confidence Intervals##
d.m.clarify<-matrix(NA,nrow=1000,ncol=12)
d.c.clarify<-matrix(NA,nrow=1000,ncol=12)
d.p.clarify<-matrix(NA,nrow=1000,ncol=12)

#Clarify simulations
drugs.media.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$m_drugs$coefficients,Sigma=vcov(drugs.VAR$varresult$m_drugs),empirical=TRUE)
drugs.cong.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$drugsCONG$coefficients,Sigma=vcov(drugs.VAR$varresult$drugsCONG),empirical=TRUE)
drugs.pres.sim<-mvrnorm(n=1000,mu=drugs.VAR$varresult$drugsPRES$coefficients,Sigma=vcov(drugs.VAR$varresult$drugsPRES),empirical=TRUE)

#Loop for simulated forecasts. This calls on standing terms defined in the mean forecast: Start endogenous variables, x00, x0, and x1.
for(j in 1:1000){
	drugs.media.coef<-drugs.media.sim[j,]
	drugs.cong.coef<-drugs.cong.sim[j,]
	drugs.pres.coef<-drugs.pres.sim[j,]

	d.m.clarify[j,1]<-m.00<-x.00%*%drugs.media.coef
	d.c.clarify[j,1]<-c.00<-x.00%*%drugs.cong.coef
	d.p.clarify[j,1]<-p.00<-x.00%*%drugs.pres.coef

	d.m.clarify[j,2]<-m.0<-x.0%*%drugs.media.coef
	d.c.clarify[j,2]<-c.0<-x.0%*%drugs.cong.coef
	d.p.clarify[j,2]<-p.0<-x.0%*%drugs.pres.coef

	all.inp<-x.1
	d.m.clarify[j,3]<-m.1<-x.1%*%drugs.media.coef
	d.c.clarify[j,3]<-c.1<-x.1%*%drugs.cong.coef
	d.p.clarify[j,3]<-p.1<-x.1%*%drugs.pres.coef
	new.pred<-c(m.1,c.1,p.1)

	for(i in 2:10){
		all.inp<-assign(paste("x",i,sep="."),c(new.pred,all.inp[-c(34:53)],1,173+i,55.398,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
		d.m.clarify[j,i+2]<-assign(paste("m",i,sep="."),get(paste("x",i,sep="."))%*%drugs.media.coef)
		d.c.clarify[j,i+2]<-assign(paste("c",i,sep="."),get(paste("x",i,sep="."))%*%drugs.cong.coef)
		d.p.clarify[j,i+2]<-assign(paste("p",i,sep="."),get(paste("x",i,sep="."))%*%drugs.pres.coef)
		new.pred<-c(get(paste("m",i,sep=".")),get(paste("c",i,sep=".")),get(paste("p",i,sep=".")))
		}
}

#Lower and Upper Bounds
d.m.5<-apply(d.m.clarify,2, quantile, .05)
d.c.5<-apply(d.c.clarify,2, quantile, .05)
d.p.5<-apply(d.p.clarify,2, quantile, .05)
d.m.95<-apply(d.m.clarify,2, quantile, .95)
d.c.95<-apply(d.c.clarify,2, quantile, .95)
d.p.95<-apply(d.p.clarify,2, quantile, .95)

#Draw Figures with Confidence intervals
Time<-c(-1:10)

##pdf('drugsMediaForecastC.pdf')
plot(y=d.m.forecast,x=Time,ylab="Media Forecast",type='l',lwd=2,ylim=c(min(d.m.5),max(d.m.95)),axes=F)
lines(y=d.m.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.m.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('May 95','July 95','Sept 95','Nov 95','Jan 95','Mar 96'))
box()
##dev.off()

##pdf('drugsCongForecastC.pdf')
plot(y=d.c.forecast,x=Time,ylab="Congressional Forecast",type='l',lwd=2,ylim=c(min(d.c.5),max(d.c.95)),axes=F)
lines(y=d.c.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.c.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('May 95','July 95','Sept 95','Nov 95','Jan 95','Mar 96'))
box()
##dev.off()

##pdf('drugsPresForecastC.pdf')
plot(y=d.p.forecast,x=Time,ylab="Presidential Forecast",type='l',lwd=2,ylim=c(min(d.p.5),max(d.p.95)),axes=F)
lines(y=d.p.5,x=Time,lty=2,col='blue',lwd=2)
lines(y=d.p.95,x=Time,lty=2,col='blue',lwd=2)
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('May 95','July 95','Sept 95','Nov 95','Jan 95','Mar 96'))
box()
##dev.off()


###robustness check with statutory indicator###
comprehensive$statutory<-comprehensive$sc_leon_drugs+comprehensive$sc_austin_drugs+comprehensive$sc_benitez_drugs+comprehensive$sc_raich_drugs
comprehensive$constitutional<-comprehensive$sc_raab_drugs+comprehensive$sc_skinner_drugs+comprehensive$sc_smith_drugs+comprehensive$sc_vernonia_drugs+comprehensive$sc_ferguson_drugs+comprehensive$sc_kyllo_drugs+comprehensive$sc_earls_drugs+comprehensive$sc_georgia_drugs+comprehensive$sc_hudson_drugs+comprehensive$sc_diaz_drugs
table(comprehensive$statutory,comprehensive$constitutional)

drugs.statutory<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$statutory)); drugs.statutory; AIC(drugs.statutory)
summary(drugs.statutory$varresult$m_drugs) 
summary(drugs.statutory$varresult$drugsCONG)
summary(drugs.statutory$varresult$drugsPRES)

drugs.constitutional<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$constitutional)); drugs.constitutional; AIC(drugs.constitutional)
summary(drugs.constitutional$varresult$m_drugs)
summary(drugs.constitutional$varresult$drugsCONG)
summary(drugs.constitutional$varresult$drugsPRES)

drugs.combined<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$statutory,comprehensive$constitutional)); drugs.combined; AIC(drugs.combined)
summary(drugs.combined$varresult$m_drugs) 
summary(drugs.combined$varresult$drugsCONG)
summary(drugs.combined$varresult$drugsPRES)

drugs.statutory.2<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$sc_leon_drugs,comprehensive$sc_austin_drugs,comprehensive$sc_benitez_drugs,comprehensive$sc_raich_drugs)); drugs.statutory.2; AIC(drugs.statutory.2)
summary(drugs.statutory.2$varresult$m_drugs)
summary(drugs.statutory.2$varresult$drugsCONG)
summary(drugs.statutory.2$varresult$drugsPRES)



###ALTERNATIVE SPECIFICATIONS SUGGESTED BY REVIEWERS###
#load alternate Mood data, merge in drug Mood
moods<-read.csv("policy_moods2015_1.csv")
drug.mood<-subset(moods, subset=subtopic==1203&year>=1981&year<=2010, select=c(year,CAPpercent))
colnames(drug.mood)[2]<-"drugMood"
comprehensive<-merge(x=comprehensive,y=drug.mood,by='year')

#load most important problem data, merge in health and crime
gallup<-read.csv("Gallups_Most_Important_Problem_3.csv")
gallup$percent<-gallup$percent*100
health.gallup<-subset(gallup, subset=majortopic==3&year>=1981&year<=2010, select=c(year,percent))
crime.gallup<-subset(gallup, subset=majortopic==12&year>=1981&year<=2010, select=c(year,percent))
colnames(health.gallup)[2]<-"healthGallup"
colnames(crime.gallup)[2]<-"crimeGallup"
two.gallup<-merge(x=health.gallup,y=crime.gallup,by="year")
comprehensive<-merge(x=comprehensive,y=two.gallup,by='year')

###Drugs with alternative version of congressional action###
#create subsets
comprehensive$no_drug_law<-comprehensive$c_bills_drugs+comprehensive$c_hearings_drugs
drugs.Alt<-subset(comprehensive, select=c(m_drugs,no_drug_law,drugsPRES))

#estimate VAR 
drugs.VAR.Alt<-VAR(drugs.Alt, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$sc_leon_drugs,comprehensive$sc_raab_drugs,comprehensive$sc_skinner_drugs,comprehensive$sc_smith_drugs,comprehensive$sc_austin_drugs,comprehensive$sc_vernonia_drugs,comprehensive$sc_ferguson_drugs,comprehensive$sc_kyllo_drugs,comprehensive$sc_earls_drugs,comprehensive$sc_benitez_drugs,comprehensive$sc_raich_drugs,comprehensive$sc_georgia_drugs,comprehensive$sc_hudson_drugs,comprehensive$sc_diaz_drugs)); drugs.VAR.Alt; AIC(drugs.VAR.Alt)

#specific tests
#cases on outcomes
summary(drugs.VAR.Alt$varresult$m_drugs) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
summary(drugs.VAR.Alt$varresult$no_drug_law) #Employment Division v. Smith (April 1990).
summary(drugs.VAR.Alt$varresult$drugsPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).
xtable(drugs.VAR.Alt$varresult$m_drugs) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
xtable(drugs.VAR.Alt$varresult$no_drug_law) #Employment Division v. Smith (April 1990).
xtable(drugs.VAR.Alt$varresult$drugsPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).

#Items on Press Coverage
a<-linearHypothesis(drugs.VAR.Alt$varresult$m_drugs,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
b<-linearHypothesis(drugs.VAR.Alt$varresult$m_drugs,c("no_drug_law.l1=0","no_drug_law.l2=0","no_drug_law.l3=0","no_drug_law.l4=0","no_drug_law.l5=0","no_drug_law.l6=0","no_drug_law.l7=0","no_drug_law.l8=0","no_drug_law.l9=0","no_drug_law.l10=0","no_drug_law.l11=0","no_drug_law.l12=0"))
c<-linearHypothesis(drugs.VAR.Alt$varresult$m_drugs,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(drugs.VAR.Alt$varresult$no_drug_law,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
e<-linearHypothesis(drugs.VAR.Alt$varresult$no_drug_law,c("no_drug_law.l1=0","no_drug_law.l2=0","no_drug_law.l3=0","no_drug_law.l4=0","no_drug_law.l5=0","no_drug_law.l6=0","no_drug_law.l7=0","no_drug_law.l8=0","no_drug_law.l9=0","no_drug_law.l10=0","no_drug_law.l11=0","no_drug_law.l12=0"))
f<-linearHypothesis(drugs.VAR.Alt$varresult$no_drug_law,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on President
g<-linearHypothesis(drugs.VAR.Alt$varresult$drugsPRES,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
h<-linearHypothesis(drugs.VAR.Alt$varresult$drugsPRES,c("no_drug_law.l1=0","no_drug_law.l2=0","no_drug_law.l3=0","no_drug_law.l4=0","no_drug_law.l5=0","no_drug_law.l6=0","no_drug_law.l7=0","no_drug_law.l8=0","no_drug_law.l9=0","no_drug_law.l10=0","no_drug_law.l11=0","no_drug_law.l12=0"))
i<-linearHypothesis(drugs.VAR.Alt$varresult$drugsPRES,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Table of Granger tests
drugs.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(drugs.granger)<-colnames(drugs.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(drugs.granger,caption="\\textbf{Granger Causality Tests for Drugs VAR Model}",label="drugs.granger",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i


###Tobacco with alternative version of congressional action###
#create subsets
comprehensive$no_tobacco_law<-comprehensive$c_bills_tobacco+comprehensive$c_hearings_tobacco
tobacco.Alt<-subset(comprehensive, select=c(m_tobacco,no_tobacco_law,tobaccoPRES))

#estimate VAR 
tobacco.VAR.Alt<-VAR(tobacco.Alt, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$sc_fda_tobacco,comprehensive$sc_lorillard_tobacco)); tobacco.VAR.Alt; AIC(tobacco.VAR.Alt)

#specific tests
#cases on outcomes
summary(tobacco.VAR.Alt$varresult$m_tobacco) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
summary(tobacco.VAR.Alt$varresult$no_tobacco_law) #Employment Division v. Smith (April 1990).
summary(tobacco.VAR.Alt$varresult$tobaccoPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).
xtable(tobacco.VAR.Alt$varresult$m_tobacco) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
xtable(tobacco.VAR.Alt$varresult$no_tobacco_law) #Employment Division v. Smith (April 1990).
xtable(tobacco.VAR.Alt$varresult$tobaccoPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).

#Items on Press Coverage
a<-linearHypothesis(tobacco.VAR.Alt$varresult$m_tobacco,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
b<-linearHypothesis(tobacco.VAR.Alt$varresult$m_tobacco,c("no_tobacco_law.l1=0","no_tobacco_law.l2=0","no_tobacco_law.l3=0","no_tobacco_law.l4=0","no_tobacco_law.l5=0","no_tobacco_law.l6=0","no_tobacco_law.l7=0","no_tobacco_law.l8=0","no_tobacco_law.l9=0","no_tobacco_law.l10=0","no_tobacco_law.l11=0","no_tobacco_law.l12=0"))
c<-linearHypothesis(tobacco.VAR.Alt$varresult$m_tobacco,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(tobacco.VAR.Alt$varresult$no_tobacco_law,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
e<-linearHypothesis(tobacco.VAR.Alt$varresult$no_tobacco_law,c("no_tobacco_law.l1=0","no_tobacco_law.l2=0","no_tobacco_law.l3=0","no_tobacco_law.l4=0","no_tobacco_law.l5=0","no_tobacco_law.l6=0","no_tobacco_law.l7=0","no_tobacco_law.l8=0","no_tobacco_law.l9=0","no_tobacco_law.l10=0","no_tobacco_law.l11=0","no_tobacco_law.l12=0"))
f<-linearHypothesis(tobacco.VAR.Alt$varresult$no_tobacco_law,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Items on President
g<-linearHypothesis(tobacco.VAR.Alt$varresult$tobaccoPRES,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
h<-linearHypothesis(tobacco.VAR.Alt$varresult$tobaccoPRES,c("no_tobacco_law.l1=0","no_tobacco_law.l2=0","no_tobacco_law.l3=0","no_tobacco_law.l4=0","no_tobacco_law.l5=0","no_tobacco_law.l6=0","no_tobacco_law.l7=0","no_tobacco_law.l8=0","no_tobacco_law.l9=0","no_tobacco_law.l10=0","no_tobacco_law.l11=0","no_tobacco_law.l12=0"))
i<-linearHypothesis(tobacco.VAR.Alt$varresult$tobaccoPRES,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Table of Granger tests
tobacco.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(tobacco.granger)<-colnames(tobacco.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(tobacco.granger,caption="\\textbf{Granger Causality Tests for tobacco VAR Model}",label="tobacco.granger",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i


###Drugs with drug policy-specific Mood##
drugs.specific.VAR<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$drugMood,comprehensive$sc_leon_drugs,comprehensive$sc_raab_drugs,comprehensive$sc_skinner_drugs,comprehensive$sc_smith_drugs,comprehensive$sc_austin_drugs,comprehensive$sc_vernonia_drugs,comprehensive$sc_ferguson_drugs,comprehensive$sc_kyllo_drugs,comprehensive$sc_earls_drugs,comprehensive$sc_benitez_drugs,comprehensive$sc_raich_drugs,comprehensive$sc_georgia_drugs,comprehensive$sc_hudson_drugs,comprehensive$sc_diaz_drugs)); drugs.specific.VAR; AIC(drugs.specific.VAR)

#cases on outcomes
summary(drugs.specific.VAR$varresult$m_drugs) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
summary(drugs.specific.VAR$varresult$drugsCONG) #Employment Division v. Smith (April 1990).
summary(drugs.specific.VAR$varresult$drugsPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).
xtable(drugs.specific.VAR$varresult$m_drugs) 
xtable(drugs.specific.VAR$varresult$drugsCONG) 
xtable(drugs.specific.VAR$varresult$drugsPRES)

#Items on Press Coverage
a<-linearHypothesis(drugs.specific.VAR$varresult$m_drugs,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
b<-linearHypothesis(drugs.specific.VAR$varresult$m_drugs,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
c<-linearHypothesis(drugs.specific.VAR$varresult$m_drugs,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(drugs.specific.VAR$varresult$drugsCONG,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
e<-linearHypothesis(drugs.specific.VAR$varresult$drugsCONG,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
f<-linearHypothesis(drugs.specific.VAR$varresult$drugsCONG,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on President
g<-linearHypothesis(drugs.specific.VAR$varresult$drugsPRES,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
h<-linearHypothesis(drugs.specific.VAR$varresult$drugsPRES,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
i<-linearHypothesis(drugs.specific.VAR$varresult$drugsPRES,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Table of Granger tests
drugs.specific.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(drugs.specific.granger)<-colnames(drugs.specific.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(drugs.specific.granger,caption="\\textbf{Granger Causality Tests for Drugs VAR Model with Policy-Specific Mood}",label="drugs.specific.granger",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i


###Drugs with crime as most important problem###
drugs.crime.VAR<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$crimeGallup,comprehensive$sc_leon_drugs,comprehensive$sc_raab_drugs,comprehensive$sc_skinner_drugs,comprehensive$sc_smith_drugs,comprehensive$sc_austin_drugs,comprehensive$sc_vernonia_drugs,comprehensive$sc_ferguson_drugs,comprehensive$sc_kyllo_drugs,comprehensive$sc_earls_drugs,comprehensive$sc_benitez_drugs,comprehensive$sc_raich_drugs,comprehensive$sc_georgia_drugs,comprehensive$sc_hudson_drugs,comprehensive$sc_diaz_drugs)); drugs.crime.VAR; AIC(drugs.crime.VAR)

#cases on outcomes
summary(drugs.crime.VAR$varresult$m_drugs) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
summary(drugs.crime.VAR$varresult$drugsCONG) #Employment Division v. Smith (April 1990).
summary(drugs.crime.VAR$varresult$drugsPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).
xtable(drugs.crime.VAR$varresult$m_drugs)
xtable(drugs.crime.VAR$varresult$drugsCONG) 
xtable(drugs.crime.VAR$varresult$drugsPRES) 


#Items on Press Coverage
a<-linearHypothesis(drugs.crime.VAR$varresult$m_drugs,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
b<-linearHypothesis(drugs.crime.VAR$varresult$m_drugs,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
c<-linearHypothesis(drugs.crime.VAR$varresult$m_drugs,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(drugs.crime.VAR$varresult$drugsCONG,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
e<-linearHypothesis(drugs.crime.VAR$varresult$drugsCONG,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
f<-linearHypothesis(drugs.crime.VAR$varresult$drugsCONG,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on President
g<-linearHypothesis(drugs.crime.VAR$varresult$drugsPRES,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
h<-linearHypothesis(drugs.crime.VAR$varresult$drugsPRES,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
i<-linearHypothesis(drugs.crime.VAR$varresult$drugsPRES,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Table of Granger tests
drugs.crime.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(drugs.crime.granger)<-colnames(drugs.crime.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(drugs.crime.granger,caption="\\textbf{Granger Causality Tests for Drugs VAR Model with Crime as Most Important Problem}",label="drugs.crime.granger",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i


###Drugs with health as most important problem###
drugs.health.VAR<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$healthGallup,comprehensive$sc_leon_drugs,comprehensive$sc_raab_drugs,comprehensive$sc_skinner_drugs,comprehensive$sc_smith_drugs,comprehensive$sc_austin_drugs,comprehensive$sc_vernonia_drugs,comprehensive$sc_ferguson_drugs,comprehensive$sc_kyllo_drugs,comprehensive$sc_earls_drugs,comprehensive$sc_benitez_drugs,comprehensive$sc_raich_drugs,comprehensive$sc_georgia_drugs,comprehensive$sc_hudson_drugs,comprehensive$sc_diaz_drugs)); drugs.health.VAR; AIC(drugs.health.VAR)

#cases on outcomes
summary(drugs.health.VAR$varresult$m_drugs) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
summary(drugs.health.VAR$varresult$drugsCONG) #Employment Division v. Smith (April 1990).
summary(drugs.health.VAR$varresult$drugsPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).
xtable(drugs.health.VAR$varresult$m_drugs) 
xtable(drugs.health.VAR$varresult$drugsCONG)
xtable(drugs.health.VAR$varresult$drugsPRES)

#Items on Press Coverage
a<-linearHypothesis(drugs.health.VAR$varresult$m_drugs,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
b<-linearHypothesis(drugs.health.VAR$varresult$m_drugs,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
c<-linearHypothesis(drugs.health.VAR$varresult$m_drugs,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(drugs.health.VAR$varresult$drugsCONG,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
e<-linearHypothesis(drugs.health.VAR$varresult$drugsCONG,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
f<-linearHypothesis(drugs.health.VAR$varresult$drugsCONG,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on President
g<-linearHypothesis(drugs.health.VAR$varresult$drugsPRES,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
h<-linearHypothesis(drugs.health.VAR$varresult$drugsPRES,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
i<-linearHypothesis(drugs.health.VAR$varresult$drugsPRES,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Table of Granger tests
drugs.health.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(drugs.health.granger)<-colnames(drugs.health.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(drugs.health.granger,caption="\\textbf{Granger Causality Tests for Drugs VAR Model with Health as Most Important Problem}",label="drugs.health.granger",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i


###Tobacco with health as most important problem###
tobacco.health.VAR<-VAR(tobacco, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$healthGallup,comprehensive$sc_fda_tobacco,comprehensive$sc_lorillard_tobacco)); tobacco.health.VAR; AIC(tobacco.health.VAR)

#cases on outcomes
summary(tobacco.health.VAR$varresult$m_tobacco) #Significant: National Treasury Employees Union v. Von Raab (March 1989) and Georgia v. Randolph (2006)
summary(tobacco.health.VAR$varresult$tobaccoCONG) #Employment Division v. Smith (April 1990).
summary(tobacco.health.VAR$varresult$tobaccoPRES) #Mood, Vernonia School District v. Acton (June 1995), and Board of Education v. Earls (June 2002).
xtable(tobacco.health.VAR$varresult$m_tobacco) 
xtable(tobacco.health.VAR$varresult$tobaccoCONG)
xtable(tobacco.health.VAR$varresult$tobaccoPRES)

#Items on Press Coverage
a<-linearHypothesis(tobacco.health.VAR$varresult$m_tobacco,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
b<-linearHypothesis(tobacco.health.VAR$varresult$m_tobacco,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
c<-linearHypothesis(tobacco.health.VAR$varresult$m_tobacco,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(tobacco.health.VAR$varresult$tobaccoCONG,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
e<-linearHypothesis(tobacco.health.VAR$varresult$tobaccoCONG,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
f<-linearHypothesis(tobacco.health.VAR$varresult$tobaccoCONG,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Items on President
g<-linearHypothesis(tobacco.health.VAR$varresult$tobaccoPRES,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
h<-linearHypothesis(tobacco.health.VAR$varresult$tobaccoPRES,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
i<-linearHypothesis(tobacco.health.VAR$varresult$tobaccoPRES,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Table of Granger tests
tobacco.health.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(tobacco.health.granger)<-colnames(tobacco.health.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(tobacco.health.granger,caption="\\textbf{Granger Causality Tests for tobacco VAR Model with Policy-health Mood}",label="tobacco.health.granger",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i



####CONTINUOUS COURT ATTENTION MEASURE####
#File obtained from http://www.comparativeagendas.net/datasets_codebooks on 30 July 2016.
court.cases<-read.csv("supreme_court_cases_2.csv")
head(court.cases)
table(court.cases$year)
#2009: 92 cases total, but only 40 in the dataset. 2010: 85 cases total, but none in the dataset. Both years omitted.

##TOBACCO CHECK##
court.cases$sc.tob<-as.numeric(court.cases$pap_subtopic==341 | court.cases$subtopic==341)
table(court.cases$sc.tob) #1945-2009, 8955 cases total, only 5 related to tobacco. They were: Cipollone v. Liggett Group (1992), FDA v. Brown & Williamson (1999), Lorillard Tobacco Company v. Reilly (2001), Rowe v. New Hampshire Motor Transport Association (2007), and Altria Group, Inc. v. Good (2008).
court.cases$case_name[court.cases$sc.tob==1]
court.cases$year[court.cases$sc.tob==1]

#Add other cases not previously considered, reestimate model
comprehensive$cipollone<-as.numeric(comprehensive$year==1992 & comprehensive$month==6)
comprehensive$rowe<-as.numeric(comprehensive$year==2008 & comprehensive$month==2)
comprehensive$altria<-as.numeric(comprehensive$year==2008 & comprehensive$month==12)
tobacco.VAR.5case<-VAR(tobacco, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$cipollone,comprehensive$sc_fda_tobacco,comprehensive$sc_lorillard_tobacco,comprehensive$rowe,comprehensive$altria)); tobacco.VAR.5case; AIC(tobacco.VAR.5case)

#specific tests
#cases on outcomes
xtable(tobacco.VAR.5case$varresult$m_tobacco) 
xtable(tobacco.VAR.5case$varresult$tobaccoCONG) 
xtable(tobacco.VAR.5case$varresult$tobaccoPRES) 

#Items on Press Coverage
a<-linearHypothesis(tobacco.VAR.5case$varresult$m_tobacco,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
b<-linearHypothesis(tobacco.VAR.5case$varresult$m_tobacco,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
c<-linearHypothesis(tobacco.VAR.5case$varresult$m_tobacco,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(tobacco.VAR.5case$varresult$tobaccoCONG,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
e<-linearHypothesis(tobacco.VAR.5case$varresult$tobaccoCONG,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
f<-linearHypothesis(tobacco.VAR.5case$varresult$tobaccoCONG,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Items on President
g<-linearHypothesis(tobacco.VAR.5case$varresult$tobaccoPRES,c("m_tobacco.l1=0","m_tobacco.l2=0","m_tobacco.l3=0","m_tobacco.l4=0","m_tobacco.l5=0","m_tobacco.l6=0","m_tobacco.l7=0","m_tobacco.l8=0","m_tobacco.l9=0","m_tobacco.l10=0","m_tobacco.l11=0","m_tobacco.l12=0"))
h<-linearHypothesis(tobacco.VAR.5case$varresult$tobaccoPRES,c("tobaccoCONG.l1=0","tobaccoCONG.l2=0","tobaccoCONG.l3=0","tobaccoCONG.l4=0","tobaccoCONG.l5=0","tobaccoCONG.l6=0","tobaccoCONG.l7=0","tobaccoCONG.l8=0","tobaccoCONG.l9=0","tobaccoCONG.l10=0","tobaccoCONG.l11=0","tobaccoCONG.l12=0"))
i<-linearHypothesis(tobacco.VAR.5case$varresult$tobaccoPRES,c("tobaccoPRES.l1=0","tobaccoPRES.l2=0","tobaccoPRES.l3=0","tobaccoPRES.l4=0","tobaccoPRES.l5=0","tobaccoPRES.l6=0","tobaccoPRES.l7=0","tobaccoPRES.l8=0","tobaccoPRES.l9=0","tobaccoPRES.l10=0","tobaccoPRES.l11=0","tobaccoPRES.l12=0"))

#Table of Granger tests
tobacco.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(tobacco.granger)<-colnames(tobacco.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(tobacco.granger,caption="\\textbf{Granger Causality Tests for Tobacco VAR Model with All Tobacco Cases}",label="tobacco.granger.5",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i


##DRUG CHECK 1##
court.cases$sc.drug<-as.numeric(court.cases$pap_subtopic==1203 | court.cases$subtopic==1203)
court.cases.2<-subset(court.cases,subset=year>=1981)
table(court.cases.2$sc.drug)
court.cases.2$case_name[court.cases.2$sc.drug==1]
court.cases.2$year[court.cases.2$sc.drug==1]

#Add other cases not previously considered, reestimate model
comprehensive$washington<-as.numeric(comprehensive$year==1982 & comprehensive$month==1)
comprehensive$hoffman<-as.numeric(comprehensive$year==1982 & comprehensive$month==3)
comprehensive$powell<-as.numeric(comprehensive$year==1984 & comprehensive$month==12)
comprehensive$hernandez<-as.numeric(comprehensive$year==1985 & comprehensive$month==7)
comprehensive$sokolow<-as.numeric(comprehensive$year==1989 & comprehensive$month==4)
comprehensive$monsanto<-as.numeric(comprehensive$year==1989 & comprehensive$month==6)
comprehensive$shabani<-as.numeric(comprehensive$year==1994 & comprehensive$month==11)
comprehensive$neal<-as.numeric(comprehensive$year==1996 & comprehensive$month==1)
comprehensive$mitchell<-as.numeric(comprehensive$year==1999 & comprehensive$month==4)
comprehensive$oakland<-as.numeric(comprehensive$year==2001 & comprehensive$month==5)
comprehensive$salinas<-as.numeric(comprehensive$year==2006 & comprehensive$month==4)
comprehensive$claiborne<-as.numeric(comprehensive$year==2007 & comprehensive$month==6)
comprehensive$gall<-as.numeric(comprehensive$year==2007 & comprehensive$month==12)
comprehensive$burgess<-as.numeric(comprehensive$year==2008 & comprehensive$month==4)
comprehensive$moore<-as.numeric(comprehensive$year==2008 & comprehensive$month==10)
comprehensive$spears<-as.numeric(comprehensive$year==2009 & comprehensive$month==1)
comprehensive$azGant<-as.numeric(comprehensive$year==2009 & comprehensive$month==4)
comprehensive$abuelhawa<-as.numeric(comprehensive$year==2009 & comprehensive$month==5)

drugs.VAR.allCase<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$washington,comprehensive$hoffman,comprehensive$sc_leon_drugs,comprehensive$powell,comprehensive$hernandez,comprehensive$sc_raab_drugs,comprehensive$sokolow,comprehensive$monsanto,comprehensive$sc_skinner_drugs,comprehensive$sc_smith_drugs,comprehensive$sc_austin_drugs,comprehensive$shabani,comprehensive$sc_vernonia_drugs,comprehensive$neal,comprehensive$mitchell,comprehensive$sc_ferguson_drugs,comprehensive$oakland,comprehensive$sc_kyllo_drugs,comprehensive$sc_earls_drugs,comprehensive$sc_benitez_drugs,comprehensive$sc_raich_drugs,comprehensive$sc_georgia_drugs,comprehensive$salinas,comprehensive$sc_hudson_drugs,comprehensive$claiborne,comprehensive$gall,comprehensive$burgess,comprehensive$moore,comprehensive$spears,comprehensive$azGant,comprehensive$abuelhawa,comprehensive$sc_diaz_drugs)); drugs.VAR.allCase; AIC(drugs.VAR.allCase)


drugs.VAR.moreCase<-VAR(drugs, ic="AIC", p=12, type="both",exogen=cbind(comprehensive$mood,comprehensive$hoffman,comprehensive$sc_leon_drugs,comprehensive$sc_raab_drugs,comprehensive$sc_skinner_drugs,comprehensive$sc_smith_drugs,comprehensive$sc_austin_drugs,comprehensive$sc_vernonia_drugs,comprehensive$neal,comprehensive$sc_ferguson_drugs,comprehensive$sc_kyllo_drugs,comprehensive$sc_earls_drugs,comprehensive$sc_benitez_drugs,comprehensive$sc_raich_drugs,comprehensive$sc_georgia_drugs,comprehensive$sc_hudson_drugs,comprehensive$gall,comprehensive$burgess,comprehensive$spears,comprehensive$sc_diaz_drugs)); drugs.VAR.moreCase; AIC(drugs.VAR.moreCase)

#specific tests
#cases on outcomes
xtable(drugs.VAR.moreCase$varresult$m_drugs) 
xtable(drugs.VAR.moreCase$varresult$drugsCONG) 
xtable(drugs.VAR.moreCase$varresult$drugsPRES) 

#Items on Press Coverage
a<-linearHypothesis(drugs.VAR.moreCase$varresult$m_drugs,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
b<-linearHypothesis(drugs.VAR.moreCase$varresult$m_drugs,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
c<-linearHypothesis(drugs.VAR.moreCase$varresult$m_drugs,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on Congress
d<-linearHypothesis(drugs.VAR.moreCase$varresult$drugsCONG,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
e<-linearHypothesis(drugs.VAR.moreCase$varresult$drugsCONG,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
f<-linearHypothesis(drugs.VAR.moreCase$varresult$drugsCONG,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Items on President
g<-linearHypothesis(drugs.VAR.moreCase$varresult$drugsPRES,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
h<-linearHypothesis(drugs.VAR.moreCase$varresult$drugsPRES,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
i<-linearHypothesis(drugs.VAR.moreCase$varresult$drugsPRES,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))

#Table of Granger tests
drugs.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2]),nrow=3,ncol=3)
rownames(drugs.granger)<-colnames(drugs.granger)<-c("Press coverage","Congressional activities","Presidential activities")
xtable(drugs.granger,caption="\\textbf{Granger Causality Tests for Drugs VAR Model}",label="drugs.granger.more",align="lrrr",digits=4)
#a;b;c;d;e;f;g;h;i






###DRUG CHECK 2###
trunc.court<-subset(court.cases,subset=year>=1981) #3469 cases.
table(trunc.court$sc.drug) #25 out of 3444 cases.
court.cases$sc.drug<-as.numeric(court.cases$pap_subtopic==1203 | court.cases$subtopic==1203)
table(court.cases$sc.drug) #1945-2009, 8955 cases total, 72 related to drugs. From 1981-2009, 0 cases were heard in 14 of the years. 
endog.court<-as.data.frame(cbind(as.vector(by(court.cases$sc.drug[court.cases$year>=1981 & court.cases$year<=2008],INDICES=court.cases$year[court.cases$year>=1981 & court.cases$year<=2008],FUN=mean)*100),c(1981:2008)))
names(endog.court)<-c('SCdrug','year')
comprehensive.b<-merge(x=comprehensive,y=endog.court,by='year')
comprehensive.b<-subset(comprehensive.b,subset=year<=2008)
drugs.b<-subset(comprehensive.b, select=c(m_drugs,drugsCONG,drugsPRES,SCdrug,mood))
drugs.VAR.b<-VAR(drugs.b, ic="AIC", p=12, type="both"); drugs.VAR.b; AIC(drugs.VAR.b)

#Items on Press Coverage
a<-linearHypothesis(drugs.VAR.b$varresult$m_drugs,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
b<-linearHypothesis(drugs.VAR.b$varresult$m_drugs,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
c<-linearHypothesis(drugs.VAR.b$varresult$m_drugs,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))
d<-linearHypothesis(drugs.VAR.b$varresult$m_drugs,c("SCdrug.l1=0","SCdrug.l2=0","SCdrug.l3=0","SCdrug.l4=0","SCdrug.l5=0","SCdrug.l6=0","SCdrug.l7=0","SCdrug.l8=0","SCdrug.l9=0","SCdrug.l10=0","SCdrug.l11=0","SCdrug.l12=0"))
e<-linearHypothesis(drugs.VAR.b$varresult$m_drugs,c("mood.l1=0","mood.l2=0","mood.l3=0","mood.l4=0","mood.l5=0","mood.l6=0","mood.l7=0","mood.l8=0","mood.l9=0","mood.l10=0","mood.l11=0","mood.l12=0"))

#Items on Congress
f<-linearHypothesis(drugs.VAR.b$varresult$drugsCONG,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
g<-linearHypothesis(drugs.VAR.b$varresult$drugsCONG,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
h<-linearHypothesis(drugs.VAR.b$varresult$drugsCONG,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))
i<-linearHypothesis(drugs.VAR.b$varresult$drugsCONG,c("SCdrug.l1=0","SCdrug.l2=0","SCdrug.l3=0","SCdrug.l4=0","SCdrug.l5=0","SCdrug.l6=0","SCdrug.l7=0","SCdrug.l8=0","SCdrug.l9=0","SCdrug.l10=0","SCdrug.l11=0","SCdrug.l12=0"))
j<-linearHypothesis(drugs.VAR.b$varresult$drugsCONG,c("mood.l1=0","mood.l2=0","mood.l3=0","mood.l4=0","mood.l5=0","mood.l6=0","mood.l7=0","mood.l8=0","mood.l9=0","mood.l10=0","mood.l11=0","mood.l12=0"))

#Items on President
k<-linearHypothesis(drugs.VAR.b$varresult$drugsPRES,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
l<-linearHypothesis(drugs.VAR.b$varresult$drugsPRES,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
m<-linearHypothesis(drugs.VAR.b$varresult$drugsPRES,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))
n<-linearHypothesis(drugs.VAR.b$varresult$drugsPRES,c("SCdrug.l1=0","SCdrug.l2=0","SCdrug.l3=0","SCdrug.l4=0","SCdrug.l5=0","SCdrug.l6=0","SCdrug.l7=0","SCdrug.l8=0","SCdrug.l9=0","SCdrug.l10=0","SCdrug.l11=0","SCdrug.l12=0"))
o<-linearHypothesis(drugs.VAR.b$varresult$drugsPRES,c("mood.l1=0","mood.l2=0","mood.l3=0","mood.l4=0","mood.l5=0","mood.l6=0","mood.l7=0","mood.l8=0","mood.l9=0","mood.l10=0","mood.l11=0","mood.l12=0"))

#Items on Supreme Court
p<-linearHypothesis(drugs.VAR.b$varresult$SCdrug,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
q<-linearHypothesis(drugs.VAR.b$varresult$SCdrug,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
r<-linearHypothesis(drugs.VAR.b$varresult$SCdrug,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))
s<-linearHypothesis(drugs.VAR.b$varresult$SCdrug,c("SCdrug.l1=0","SCdrug.l2=0","SCdrug.l3=0","SCdrug.l4=0","SCdrug.l5=0","SCdrug.l6=0","SCdrug.l7=0","SCdrug.l8=0","SCdrug.l9=0","SCdrug.l10=0","SCdrug.l11=0","SCdrug.l12=0"))
t<-linearHypothesis(drugs.VAR.b$varresult$SCdrug,c("mood.l1=0","mood.l2=0","mood.l3=0","mood.l4=0","mood.l5=0","mood.l6=0","mood.l7=0","mood.l8=0","mood.l9=0","mood.l10=0","mood.l11=0","mood.l12=0"))

#Items on Public Opinion
u<-linearHypothesis(drugs.VAR.b$varresult$mood,c("m_drugs.l1=0","m_drugs.l2=0","m_drugs.l3=0","m_drugs.l4=0","m_drugs.l5=0","m_drugs.l6=0","m_drugs.l7=0","m_drugs.l8=0","m_drugs.l9=0","m_drugs.l10=0","m_drugs.l11=0","m_drugs.l12=0"))
v<-linearHypothesis(drugs.VAR.b$varresult$mood,c("drugsCONG.l1=0","drugsCONG.l2=0","drugsCONG.l3=0","drugsCONG.l4=0","drugsCONG.l5=0","drugsCONG.l6=0","drugsCONG.l7=0","drugsCONG.l8=0","drugsCONG.l9=0","drugsCONG.l10=0","drugsCONG.l11=0","drugsCONG.l12=0"))
w<-linearHypothesis(drugs.VAR.b$varresult$mood,c("drugsPRES.l1=0","drugsPRES.l2=0","drugsPRES.l3=0","drugsPRES.l4=0","drugsPRES.l5=0","drugsPRES.l6=0","drugsPRES.l7=0","drugsPRES.l8=0","drugsPRES.l9=0","drugsPRES.l10=0","drugsPRES.l11=0","drugsPRES.l12=0"))
x<-linearHypothesis(drugs.VAR.b$varresult$mood,c("SCdrug.l1=0","SCdrug.l2=0","SCdrug.l3=0","SCdrug.l4=0","SCdrug.l5=0","SCdrug.l6=0","SCdrug.l7=0","SCdrug.l8=0","SCdrug.l9=0","SCdrug.l10=0","SCdrug.l11=0","SCdrug.l12=0"))
y<-linearHypothesis(drugs.VAR.b$varresult$mood,c("mood.l1=0","mood.l2=0","mood.l3=0","mood.l4=0","mood.l5=0","mood.l6=0","mood.l7=0","mood.l8=0","mood.l9=0","mood.l10=0","mood.l11=0","mood.l12=0"))

#Table of Granger tests
drugs.granger<-matrix(c(a$F[2],b$F[2],c$F[2],d$F[2],e$F[2],f$F[2],g$F[2],h$F[2],i$F[2],j$F[2],k$F[2],l$F[2],m$F[2],n$F[2],o$F[2],p$F[2],q$F[2],r$F[2],s$F[2],t$F[2],u$F[2],v$F[2],w$F[2],x$F[2],y$F[2]),nrow=5,ncol=5)
rownames(drugs.granger)<-colnames(drugs.granger)<-c("Press coverage","Congressional activities","Presidential activities","Supreme Court","Mood")
xtable(drugs.granger,caption="\\textbf{Granger Causality Tests for Drugs VAR Model with the Court and Public Opinion as Endogenous}",label="drugs.granger",align="lrrrrr",digits=4)
#a;b;c;d;e;f;g;h;i;j;k;l;m;n;o;p;q;r;s;t;u;v;w;x;y





